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DISCLAIMER 


Neither  the  developer  nor  any  person  or  organization  acting  on 
behalf  of  him: 

a)  makes  any  warranty,  express  or  implied,  with  respect  to  this 
software;  or 

b)  assumes  any  liabilities  with  respect  to  the  use,  or  misuse,  of 
this  software,  or  the  interpretation,  or  misinterpretation,  of 
any  results  obtained  from  this  software,  or  for  damages  resulting 
from  the  use  of  this  software. 


ABSTRACT 


GWPATH  is  an  interactive  software  package  for  estimating 
horizontal  fluid  pathlines  and  travel  times  in  fully  saturated  ground- 
water  flow  domains.  GWPATH  was  developed  for  the  IBM  PG-AT,  and 
compatibles,  microcomputing  environment  and  takes  full  advantage  of 
interactive  graphical  display.  The  code  features  interactive  data 
entry,  forward  and  reverse  pathline  tracking,  time -related  capture  zone 
analysis,  multiple  pathline  capture  detection  mechanisms,  multiple 
configurations  for  pathline  starting  locations,  and  variable  time 
stepping.  GWPATH  is  applicable  to  inhomogeneous,  anisotropic  complex 
flow  domains.  Data  requirements  include  definition  of  the  ground-water 
flow  domain  configuration,  a  hydraulic  head  distribution,  hydraulic 
conductivity,  and  effective  porosity.  Both  hardcopy  printable  output 
and  six-pen  plotter  output  are  generated  by  GWPATH. 


TABLE  OF  NOTATION 


1.  Mathematical  Symbols 


h 

q 

K 

f(x) 

g(x,y) 

da/dx 

da/dx 

grad  h 


A  nonboldface  lower  case  letter:  a  scalar. 

A  boldface  lower  case  letter:  a  vector  in  two  dimensions  with 
components  and  qy. 

A  boldface  upper  case  letter:  a  tensor. 

A  function  of  one  variable. 

A  function  of  two  variables. 

The  first  order  derivative  of  a  =  f(x)  with  respect  to  x. 

The  first  order  partial  derivative  of  a  =  g(x,y)  with  respect 
to  X  . 

(5h/3x,  dh/dy) :  the  gradient  of  scalar  h. 


2.  Variables 


h 

K 


q 

s(t) 

V 


x,y 


Total  hydraulic  head  (L) . 

Hydraulic  conductivity  tensor  oriented  along  principal 
horizontal  axes  of  anisotropy  (L/T) . 


Effective  porosity  (dimensionless) . 

Darcy  velocity  vector  with  two-dimensional  horizontal 
components  qjr  and  qy  (L/T)  . 

Ground-water  pathline  trajectory:  s(t)  =  [x(t),y(t)]  (L) . 
Average  pore  ground-water  velocity  with  two-dimensional 
horizontal  components  v^  and  Vy  (L/T) . 

Horizontal  position  coordinates  (L) . 


PROGRAM  SYNOPSIS 


Program  Name:  GWPATH,  Version  3.0 

Purpose:  Estimate  two-dimensional  horizontal  ground-water  flow  paths  in 
complex  flow  domains.  Track  flow  paths  down- gradient  from 
source  to  sink;  or  up-gradient  from  sink  to  source.  Compute 
ground-water  travel  time  (days)  for  each  pathline.  Estimate 
time-related  capture  zones  for  hydraulic  sinks. 

Approach:  Lagrangian 

Language:  IBM  Professional  FORTRAN,  Version  1.14 

Mode:  Interactive 


Recommended  Hardware  :  IBM  AT  w/  hard  disk  /  math  co-processor 

EGA  /  high  resolution  color  monitor 
HP-7475A  or  IBM-7372  six-pen  plotter 

Software  Requirements:  DOS  2.x  or  greater 

Input  Requirements:  Flow  domain  configuration  (i.e.,  grid  structure) 

Hydraulic  head  distribution 
Hydraulic  conductivity 
Effective  porosity 
Source/sink  locations  (optional) 

Output:  Screen  graphics  display 

Disk  file  of  pathline(s)  trajectory  and  travel  time 
Optional  plot  (HP-7475A/IBM-7372  six-pen  plotter) 

Assumptions:  Regular,  nonuniform,  mesh- centered  rectangular  grid 
Steady  state  hydraulic  head  distribution 
Two-dimensional  regional  flow 

Features:  Interactive  data  entry 

Inhomogeneous,  anisotropic  flow  domain 
Continuous  velocity  distribution 
Fourth-order  Runge-Kutta  numerical  integration 
Variable  time  stepping 

Multiple  configurations  for  pathline  starting  locations 
Forward  tracking;  reverse  tracking 
Multiple  pathline  capture  detection  mechanisms 
Time-related  capture  zone  analysis 


Note:  GWPATH  supports  only  IBM  standard  devices. 
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METHODOLOGY 


The  following  sections  describe  the  methodology  implemented  in 
GWPATH  for  estimation  of  two-dimensional,  horizontal  ground-water  flow 
paths.  For  more  in-depth  discussions  of  ground-water  hydraulics  and/or 
ground-water  flow  modeling,  several  informative  general  references  are 
available  including: 

Bear,  J.,  1979.  Hydraulics  of  Groundwater,  McGraw-Hill  Book  Co., 

New  York,  New  York,  569  pp . 

de  Marsily,  G.,  1986.  Quantitative  Hydrogeology,  Academic  Press, 

Inc.,  Orlando,  Florida,  440  pp . 

Mercer,  J.W.  and  C.R.  Faust,  1981.  Ground-Water  Modeling, 

National  Water  Well  Association,  Worthington,  Ohio,  60  pp . 

Wang,  H.F.  and  M.P.  Anderson,  1982.  Introduction  to  Ground- 
water  Modeling,  W.H.  Freeman  and  Co.,  San  Francisco, 
California,  237  pp . 


Overview 

There  are  four  basic  types  of  ground-water  models,  each  addressing 
different  physical  phenomena;  i.e.,  ground-water  flow,  solute  transport, 
heat  transport,  and  structural  deformation  (Mercer  and  Faust,  1981). 
Ground-water  pathline  analysis  is  an  extension  of  ground-water  flow 
modeling.  It  seeks  a  balance  between  predictive  sophistication  (i.e., 
solute  transport  modeling)  and  the  amount  and  quality  of  available  data. 
Solute  concentrations  are  not  directly  modeled;  however,  center-of -mass 
flow  path  trajectories  and  ground-water  travel  times  are  modeled.  There 
are  a  variety  of  regional  and  sub-regional  ground-water  contamination 
problems  where  flow  path  determination  (especially  for  contaminant 
source  identification)  is  a  logical  first  step  of  analysis  (Nelson  and 
Schur,  1980).  Further,  pathline/travel  time  analyses,  used  as 
precursors  to  more  elaborate  evaluations,  can  aid  in  optimizing  data 
collection  efforts. 

The  reverse  pathline  calculation  capability  of  GWPATH  provides  the 
capability  for  estimating  time-related  capture  zones  of  hydraulic  sinks 
throughout  the  flow  domain.  A  large  number  (i.e.,  usually  between  100 
and  300)  of  reverse  calculated  pathlines,  all  with  the  same  travel  time, 
can  be  used  to  approximate  the  envelope  of  influence  of  hydraulic  sinks 
over  the  specified  time  period.  For  a  detailed  discussion  of  reverse 
pathline  calculation  of  time-related  capture  zones,  refer  to  Shafer 
(1987)  . 
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Theoretical  Background 

The  ground-water  flow  system  kinematics  provide  the  geometry  of 
saturated  ground-water  movement.  A  pathline  is  the  projected 
two-dimensional  representation  of  the  mean  trajectory  and  travel  time  an 
individual  fluid  particle*  follows  through  the  flow  domain.  The 
pathline  describes  the  expected  planar  location  and  time  of  travel  of 
the  fluid  from  any  initial  location,  along  its  trajectory,  to  its 
terminus  (e.g.,  well  or  flow  domain  boundary).  Trajectories,  or 
pathlines,  can  be  tracked  forward  (i.e.,  down- gradient )  or  in  reverse 
(i.e.,  up-gradient)  through  the  flow  system.  Multiple  pathlines  can  be 
used  to  determine  the  envelope  of  influence  of  sinks  within  the  flow 
domain;  i.e.,  the  time-related  capture  zone. 

The  basic  solution  approach  used  to  calculate  pathlines  is  La- 
grangian.  An  infinitesimal  mass  of  fluid  (the  particle)  is  followed  as 
it  would  move  through  the  horizontal,  two-dimensional  representation  of 
the  flow  domain.  The  equations  for  fluid  pathlines  are  derived  from 
Darcy's  law  for  incompressible  fluids  (see  de  Marsily,  1986),  which  is 

q  =  -K  grad  h  (1) 


where 

q  =  Darcy  velocity  vector  with  two-dimensional  horizontal 
components  q^  and  qy  (L/T) , 

K  =  hydraulic  conductivity  tensor  with  two-dimensional 

components  and  Ky,  assuming  the  coordinate  system  is 
aligned  with  the  principal  horizontal  axes  of  anisotropy 
(L/T), 

h  =  total  hydraulic  head,  h  =  h(x,y)  (L) ,  and 
x,y  =  two-dimensional  horizontal  Cartesian  coordinates  (L) . 


The  average  pore  ground-water  velocity  is 


where 


V 


=  3_ 


n 


e 


V  =  average  pore  velocity  vector  with  two-dimensional 
horizontal  components  v^,  and  Vy  (L/T)  ,  and 
ng  =  effective  porosity  (dimensionless) . 


(2) 


Combining  equations  1  and  2  results  in 


v 


grad  h 


(3) 


This  "particle"  has  no  physical  meaning  and 
mathematical  representation  for  an  infinitesimal  mass  of 


is  used 
fluid. 


as 


a 
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with  individual  horizontal  velocity  components 


f(x,y)  =  =  Vx  (4.1) 

and  , 

g(x,y)  =  ^  =  ^y  (^-2) 

A  pathline,  s(t)  =  [x(t),y(t)],  in  the  flow  domain  is  one  for 

which  the  velocity,  v,  associated  with  the  path  satisfies  equation  3 
everywhere  along  the  path.  The  determination  of  pathlines  is  achieved 
by  simultaneously  solving  the  initial  value  problems  expressed  by  the 
set  of  ordinary  differential  equations  (equations  4.1  and  4.2),  which 
represent  an  equivalent  form  of  equation  3.  Initial  conditions  are 

x(0)  =  Xq 

y(0)  =  yo- 


Numerical  Implementation 

The  solution  for  two-dimensional  ground-water  fluid  pathlines  is 
accomplished  by  numerically  implementing  the  approach  outlined  in  the 
previous  section.  There  are  three  basic  components  to  the  numerical 
determination  of  pathlines:  1)  calculation  of  Lnternodal  average  pore 
ground-water  velocities,  2)  interpolation  of  intracell  pore  ground- 
water  velocities,  and  3)  numerical  integration  of  equations  4.1  and  4.2. 
The  methods  used  for  each  component  are  discussed  in  the  following 
paragraphs.  GWPATH  is  capable  of  solving  for  pathlines  in 
inhomogeneous,  anisotropic  flow  domains.  Consequently,  two  sets  of 
internodal  average  pore  ground-water  velocities  are  calculated.  The 
first  set  is  velocities  between  nodes  that  are  horizontally  adjacent 
(i.e.,  west  to  east  in  the  x-direction)  throughout  the  flow  domain.  The 
second  set  of  velocities  is  between  nodes  that  are  vertically  adjacent 
(i.e.,  north  to  south  in  the  y-direction) .  If  the  forward  (i.e., 
down- gradient )  mode  is  selected,  all  internodal  velocities  are 
multiplied  by  -1.0.  Conversely,  if  reverse  (i.e.,  up-gradient)  tracking 
is  desired,  all  internodal  velocities  are  multiplied  by  +1.0. 

Figure  1  shows  the  mesh-centered  flow  domain  grid  configuration 
with  orientation  of  flow  parameters.  The  hydraulic  conductivities  input 
to  CWPATH  are  assumed  to  be  pointwise  directional,  measured  at  the 
nodes.  An  arithmetic  mean  is  computed  to  represent  the  internodal 
hydraulic  conductivity.  Therefore,  two  sets  of  hydraulic  conductivities 
are  used;  i.e.,  one  set  in  the  x-direction  and  one  set  in  the 
y-direction.  Effective  porosities  are  assvuned  to  be  representative  of 
the  porosity  at  each  node  and  are  scalar  quantities. 


Internodal  Velocity  Calculation.  The  method  used  to  calculate  the 
internodal  pore  ground-water  velocities  in  the  x-direction  is 
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I 


i-i 


i+i 


◄-DELX(I)-^ 

nj) 

+' 

> 

UJ 

rv 

i . 

I 

■hdj) 

Kx(lJ) 

Ky(lJ) 

ne(IJ) 

Figure  1.  Flow  domain  grid  configuration  and  orientation  of  ground- 
water  flow  parameters 


where 


VELX(i,j) 


Kx(i-.j) 

5h(i,j)' 

.ne(i.j). 

.  5x(i)  ^ 

(5) 


Kx(iJ)  =  [Kx(i+l,j)  +  Kxd.j)]  /  2 
ne(i,j)  “  (ne(i+l,j)  +  ne(i,j)]  /  2 
3h(i,j)  »  [h(i+l,j)  -  h(i,j)] 
ax(i)  =  DELX(i) 
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Similarly,  in  the  y-direction  internodal  velocities  are  calculated  as 


where 


VELY(i,j) 


Ky(i,j)^ 

'3h(i,j)' 

ne(i.j). 

.  3y(j)  ^ 

(6) 


Ky(i,j)  =  [Ky(i,j+1)  +  Ky(i,j)]  /  2 
ne(i-.j)  =  [ne(i,j+l)  +  ne(i,j)]  /  2 
ah(i,j)  =  [h(i,j+l)  -  h(i,j)] 
ay(i)  =  DELY(j) 

The  internodal  velocities  are  calculated  once  and  saved.  They  are  the 
basis  for  interpolation  of  intracell  velocities  throughout  the  entire 
flow  domain. 

Intracell  Velocity  Interpolation.  The  method  of  interpolating 
intracell  velocity  components  at  any  (x,y)  location  within  the  flow 
domain  is  based  on  an  approach  developed  by  Prickett,  Naymik,  and 
Lonnquist  (1981) .  It  is  an  elaborate  scheme  created  to  maintain 
continuity  of  velocity  across  cell  boundaries.  A  continuous  velocity 
field  is  important  in  assuring  high  accuracy,  especially  for  divergent 
flow  which  occurs,  for  example,  away  from  the  center  of  an  injection 
well.  Prickett,  Naymik,  and  Lonnquist  (1981)  developed  a  three-step 
procedure  to  interpolate  velocities  using  linear  basis  functions  (i.e.. 
Chapeau  functions)  borrowed  from  finite  element  theory. 

The  first  interpolation  uses  the  four  immediate  internodal 
velocities  surrounding  the  current  (x,y)  location.  These  are  the 
outlined  arrows  in  figure  2.  The  relative  position  within  a  grid  cell 
is  used  to  calculate  weighted  averages  of  x-direction  and  y-direction 
internodal  velocities  for  the  initial  interpolated  velocities,  VX  and 
VY,  at  (x,y). 

The  second  interpolation  uses  the  four  next  closest  internodal 
velocities  (figure  3).  For  the  demonstration  presented  in  figure  3, 
these  velocities  are  the  outlined  arrows.  In  exactly  the  same  manner  as 
the  first  interpolation,  weighted  averages,  based  on  the  (x,y)  position, 
are  used  to  calculate  the  two  intermediate  velocity  components,  vk  and 
i/y. 


Finally,  the  third  step  combines  the  first  and  second 
interpolations  (figure  4)  .  The  a  and  b  coefficients  from  the  linear 
basis  functions  are  used  to  combine  VX  and  i/x  (VY  and  uy)  to  estimate 
the  intracell  velocity  components  at  position  (x,y).  With  the  x- 
direction  for  explanation,  XZ  represents  a  weighting,  or  proportioning, 
between  VX  and  i/x .  As  the  position  x  nears  I,  XZ  increases  to  a  maximum 
of  0.5  as  b  goes  to  zero  at  I.  Hence,  the  interpolated  velocity  is 
continuous  across  I.  The  same  procedure  is  repeated  for  the  y-direction 
velocity  component.  This  process  results  in  a  continuous  velocity  field 
across  cell  boundaries.  Figure  5  shows  the  internodal  velocity 
components  used  in  the  interpolation  scheme  based  on  location  within  the 
flow  domain. 
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VELOCITY  VECTORS  BETWEEN 


DELX(I) 


COLUMNS  I  AND  I+l 


>- 


CHAPEAU  FUNCTIONS  FOR 
VELOCITY  VECTORS  BETWEEN 
ROWS  J  AND  J  +  1 


VX  =  aVELX(I,J  +  l)  + 

(l-a)VELX(IJ) 

VY  =  bVELY(I+l,J)  + 

(l-b)VELY(IJ) 


Figure  2.  Initial  interpolation  of  intracell  pore  ground-water 

velocities  (after:  Prickett,  Naymik,  and  Lonnquist,  1981) 
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VELOCITY  VECTORS  BETWEEN 
COLUMNS  I-l  AND  I 

DELX(I) 

aVELX(I-l,J  +  l) 

+  (l-a)VELX(I-l,J) 

bVELY(I+lJ-l) 

+  (l-b)VELY(IJ-l) 


I-l  I  I+l 

CHAPEAU  FUNCTIONS  FOR 
VELOCITY  VECTORS  BETWEEN 
ROWS  J-1  AND  J 


+  1 

C’ 

0 


ux  = 


uy  = 


Figure  3.  Intermediate  interpolation  of  intracell  pore  ground-water 
velocities  (after:  Prickett,  Naymik,  and  Lonnquist,  1981) 
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CHAPEAU  FUNCTIONS  FOR 
y  DIRECTION  VELOCITY  VECTORS 


XZ  1-XZ 


CHAPEAU  FUNCTIONS  FOR 
X  DIRECTION  VELOCITY  VECTORS 


XZ  =  |b-0.5| 

YZ  =  |a-0.5| 

FINAL  VELOCITY  COMPONENTS 
VX  =  (l-XZ)VX  +  XZ-ux 
VY  =  (l-YZ)VY  +  YZ'uy 


Figure  4.  Final  interpolation  of  intracell  pore  ground-water 

velocities  (after;  Prickett,  Naymik,  and  Lonnquist,  1981) 
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>•  X 

INTERIOR  OF  MODEL 


M  I  I  +  l 


J=1 


SIDES  OF  MODEL 

1=1 


I=NC 


CORNERS  OF  MODEL 


Figure  5.  Internodal  velocity  components  used  in  GWPATH 

interpolation  scheme  (after:  Prickett,  Naymik,  and 
Lonnquist,  1981) 
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Numerical  Integration.  The  determination  of  fluid  pathlines  is 
achieved  by  numerically  solving  the  linked  system  of  equations  (i.e., 
equations  4.1  and  4.2)  using  a  fourth  order  Runge-Kutta  numerical 
integration  scheme.  There  are  numerous  references  to  this  approach  in 
the  literature.  A  discussion  specific  to  pathline  analysis  can  be  found 
in  Daly  and  Morel-Seytoux  (1980).  A  variably  calculated  time  step,  At 
(days),  is  used  to  calculate  successive  locations  [x(n),y(n)], 
[x(n+l) ,y (n+1) ]  according  to 

x(n+l)  =  x(n)  +  l/6[ap  +  2a2  +  2a3  +  ai^]  (7.1) 

y(n+l)  =  y(n)  +  l/6[b3  +  2b2  +  2b3  +  b4]  (7.2) 

where 


a^  =  At  X  f [x(n) ,y(n) ] 
bx  =  At  X  g[x(n) ,y(n) ] 
a2  =  At  X  f [x(n)+ax/2 , y (n)+bx/2 ] 
b2  =  At  X  g[x(n)+ax/2 ,y(n)+bx/2] 
a3  =  At  X  f [x(n)+a2/2 , y (n)+b2/2 ] 
b3  =  At  X  g[x(n)+a2/2 ,y(n)+b2/2] 
a4  =  At  X  f [x(n)+a3 ,y(n)+b3] 
b4  =  At  X  g[x(n)+a3,y(n)+b3] 


The  Runge-Kutta  algorithm  requires  the  starting  location  (xQ,yQ). 


Time  Stepping.  The  time  step.  At,  is  allowed  to  vary  according  to 
user-defined  criteria.  At  the  beginning  of  each  Runge-Kutta  iteration 
the  time  step  is  calculated  to  be  approximately  equal  to  the  time 
required  for  the  fluid  to  linearly  traverse  the  longest  dimension  of  the 
current  cell  divided  by  the  user- supplied  number  of  moves  per  cell.  The 
larger  the  number  of  moves,  the  more  accurate  the  solution  (i.e.,  the 
shorter  the  time  step).  However,  the  shorter  the  time  step,  the  greater 
the  computational  time  necessary  to  reach  a  solution.  The  user  is  also 
required  to  enter  the  maximum  desired  time  step  which  is  invoked  if  the 
variable  time  step  exceeds  the  maximum  allowable  At.  All  travel  times 
are  in  days . 


Pathline  Capture.  Pathlines  are  often  captured  by  sinks  (forward 
tracking)  or  sources  (reverse  tracking)  within  the  flow  domain.  If 
numerical  capture  mechanisms  were  not  included  in  the  model  formulation, 
pathlines  would  oscillate  indefinitely  around  sinks  (or  sources) .  A 
dual  approach  is  used  in  GWPATH  to  terminate  pathlines  at  capture  points 
with  travel  times  accurate  to  At.  First,  point  locations  (i.e.,  x,y 
coordinates)  of  all  identifiable  sinks  and  sources  should  be  determined 
for  the  flow  domain.  These  locations  are  part  of  the  optional  input  to 
GWPATH.  At  the  beginning  of  each  time  step,  GWPATH  calculates  the 
distances  from  the  current  pathline  position  to  the  location  of  each 
sink  (or  source).  GWPATH  then  projects  the  travel  distance  of  the 
pathline  using  the  interpolated  Vj^  and  Vy  velocity  components  associated 
with  the  current  position.  If  the  projected  travel  distance  is  less 
than  the  distance  to  all  sinks  (or  sources),  the  trajectory  is  extended 
At  with  no  capture.  However,  if  the  projected  pathline  travel  distance 
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I 


I+l 


di  >  di  ;Continue 

d2  <  di  ;Backup  time  step,  reduce  A  t. 


Figure  6.  Schematic  representation  of  pathline  capture  detection 
procedure 
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is  greater  than  one  or  more  of  the  distances  to  each  sink  (or  source)  , 
the  pathline  is  assumed  to  be  captured  by  the  sink  (or  source)  closest 
to  the  current  pathline  location. 

GWPATH  has  an  override  feature  to  the  pathline  capture  detection 
method  described  above,  which  is  diagrammatically  presented  in  figure  6. 
If  during  any  particular  time  step  the  distance  from  the  previous 
position  [x(n-l) ,y(n-l) ]  to  the  next  position  [x(n+l) , y (n+1) ]  is  less 
than  the  distance  from  the  previous  position  to  the  current  position 
[x(n),y(n)]  a  flag  is  raised.  This  situation  indicates  that  there  has 
probably  been  a  90  degree  or  more  reversal  in  the  direction  of  flow. 
The  dashed  line  in  figure  6  is  orthogonal  to  the  current  trajectory. 
Such  an  occurrence  is  indicative  of  trajectory  oscillation  around  a  sink 
(or  source).  In  this  event,  the  entire  numerical  integration  procedure 
is  backed  up  one  time  step  and  a  new  mode  is  entered.  The  time  step  is 
reduced  by  10  percent  and  the  trajectory  is  recalculated  beginning  at 
location  [x(n- 1) , y (n- 1) ] .  This  process  is  repeated,  each  time  reducing 
the  time  step  by  10  percent,  until  either  the  flag  is  no  longer  raised 
or  the  time  step  has  been  reduced  to  less  then  one  minute.  If  At  has 
been  reduced  to  less  than  one  minute,  the  terminus  of  the  pathline  is 
assumed  to  be  [x(n-l) ,y(n-l) ] . 


As sump tions/Limi tat ions 

There  are  several  tacit  assumptions  associated  with 
two-dimensional  numerical  pathline  analysis.  The  ground-water  fluid 
flow  paths  produced  by  GWPATH  are  projected  average  pathlines.  In 
reality  ground  water  flows  in  three  dimensions.  With  two-dimensional 
horizontal  pathline  analysis  the  vertical  flow  component  is  assumed 
zero.  This  implies  that  along  any  vertical  section  of  the  flow  domain 
all  particles  are  moving  horizontally  in  the  same  direction  and  at  the 
same  speed. 

As  configured,  GWPATH  assumes  that  the  hydraulic  head  distribution 
is  steady  throughout  the  time  plane  of  pathline  calculation.  However, 
successive  time  planes  can  be  averaged  to  approximate  the  impact  of  a 
time-varying  hydraulic  head  distribution  on  flow  path  trajectories.  In 
this  sense,  the  resulting  pathlines  become  an  averaged  estimate  of 
trajectories  and  travel  times,  not  a  true  approximation  of  a  transient 
traj  ectory . 

As  previously  discussed,  GWPATH  requires  a  regular,  although  not 
necessarily  uniform,  mesh-centered  rectangular  grid  to  operate.  Large 
irregularly  shaped  flow  domains  can  be  subset  into  regular  components 
with  GWPATH  implemented  on  the  various  domain  subsets.  An  important 
consideration  is  that  the  accuracy  of  the  calculated  pathlines  is  a 
function  of  the  quality  and  quantity  of  available  data  as  well  as  the 
grid  density  and  length  of  time  step. 
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INPUT  REQUIREMENTS 


GWPATH  is  designed  for  interactive  data  input  via  query  response. 
The  code  is  organized  so  that  default  uniform  conditions  (i.e., 
homogeneous  hydraulic  conductivity,  effective  porosity,  and  flow  domain 
discretization)  can  be  input  directly  or  variable  conditions  can  be  read 
from  disk  files.  Hydraulic  head  distributions  are  input  from  ASCII 
files.  Figure  7  shows  the  series  of  data  input  screens. 


Units  of  measure  are  user-defined  as  either  U.S.  Customary  or  SI. 
Time  units  are  always  days.  The  following  consistencies  must  be 
observed . 


U.S.  Customary  units 
Time:  days 

Flow  domain  definition: 
Hydraulic  conductivity: 
Velocity:  feet/day 
Pathline:  feet 
SI  units 

Time:  days 

Flow  domain  definition: 
Hydraulic  conductivity: 
Velocity:  meters/day 
Pathline:  meters 


feet 

feet/day 


meters 

meters/day 


2-D  HORIZONTAL  GROUND-WATER  PATHLINE  ANALYSIS 
Version  3.0,  1987 


DISCLAIMER 

Neither  the  developer  nor  any  person  acting  on  behalf  of  him: 

(a)  makes  any  warranty,  express  or  implied,  with  respect  to  this  software;  or 

(b)  assumes  any  liabilities  with  respect  to  the  use,  or  misuse,  of  this 
software,  or  the  interpretation,  or  misinterpretation,  of  any  results 
obtained  from  this  software,  or  for  damages  resulting  from  the  use  of 
this  software. 


Execution  suspended  :  Press  ENTER  to  Resume  Execution. 


Figure  7a.  Screen  1:  Program  title  and  disclaimer 
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There  are  five  categories  of  required  data  input  for  GWPATH: 

1.  Title  and  units 

2.  Pathline  domain  definition 

3.  Flow  parameters 

4.  Output  file  specifications 

5.  Program  execution  control  information 

Ground-Water  Flow  Domain  Definition 

The  flow  domain  is  always  assumed  oriented  so  that  the  positive 
x-direction  is  from  west  to  east,  and  the  positive  y-direction  is  from 
north  to  south  on  a  flat  page.  Therefore,  the  location  of  the  minimum  x 
[XMIN]  ,  minimum  y  [YMIN]  dimension  is  in  the  upper  left  corner.  The 
user  is  asked  to  input  XMIN  and  YMIN,  the  number  of  column  nodes  [NC], 
and  the  number  of  row  nodes  [NR]  of  the  flow  domain.  The  user  is  also 
prompted  for  x-direction  spacing  [DELX]  and  y-direction  spacing  [DELY] . 
The  fields  provided  for  response  must  be  completely  filled. 

If  data  are  to  be  entered  from  files,  pressing  <Enter>  (instead  of 
typing  a  numeric  value)  will  reveal  the  filename  input  query.  In  this 
manner  a  nonuniform  grid  can  be  used  to  describe  the  ground-water  flow 


GWPATH  Version  3.0,  1987 


Title  For  This  Session  [This  is  an  Input  Demonstration] 
Units  of  Length  (Feet  or  Meters)  [Feet  ] 


PATHLINE  DOMAIN  DEFINITION 


West  (XMIN)  Boundary  [0.000000] 
No.  of  Column  Nodes  (NC)  [51] 


North  (YMIN)  Boundary  [0.000000] 
No.  of  Row  Nodes  (NR)  [51] 


GRID  SPACING 


Enter  Numeric  Real  Value  If  Uniform,  Otherwise 
Press  ENTER  For  Variable  Spacing  Input  File 


Delta-X,  East-West  Spacing  [160.00]  Delta-Y,  North-South  Spacing  [160.00] 


Figure  7b.  Screen  2:  Flow  field  definition  input 
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domain.  Variable  x-direction  and/or  y-direction  flow  domain  grid 
spacing  may  be  input  from  ASCII  files  constructed  for  list-directed 
input  according  to  the  following  example. 


Record  No. 

1 

2 

3 

I 

NC-1 

EOF 


Value 

DELX(l) 

DELK(2) 

DELX(3) 

DELX(I) 

DELX  (NC-1) 


There  must  be  NC-1  entries  in  the  DELX  file  and  NR-1  entries  in  the  DELY 
file . 


Ground-Water  Flow  Parameters 


Hydraulic  conductivity  data  and  effective  porosity  data  are 
entered  in  a  manner  identical  to  pathline  domain  grid  definition  input. 


GWPATH  Version  3.0,  1987 
FLOW  PARAMETERS 

Enter  Numeric  Real  Value  If  Uniform,  Otherwise 
Press  ENTER  For  Variable  Parameter  Input  File 


Hydraulic  Conductivity  (feet/day) 

X-Direction  [50.000]  Y-Direction  [  ] 

File  Name  [HYDY1.DAT  ] 


Effective  Porosity  (%)  [15] 


File  Containing  Hydraulic  Head  Distribution  [GWHEAD1.DAT  ] 


File  Containing  Source/Sink  Information  [SSI. DAT  ] 

Press  ENTER  If  None 


Figure  7c.  Screen  3:  Flow  parameter  input 
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If  hydraulic  conductivity  is  uniform  in  either  the  x-  or  y-direction  the 
actual  value  is  input  (the  field  must  be  completely  filled).  For 
inhomogeneous  conditions,  pressing  <Enter>  will  reveal  the  filename 
query.  The  name  of  a  disk  file  that  contains  the  variable  hydraulic 
conductivity  data  is  entered.  The  same  procedure  is  used  for  effective 
porosity.  ASCII  grid  data  matrices  must  be  formatted  using  a  ROW, 
BEGINNING  COLUMN,  ENDING  GOLUMN  arrangement  with  up  to  six  grid-node 
values  per  80  column  record.  In  this  manner,  a  typical  record  would 
have  the  matrix  row  number  for  the  first  entry  followed  by  a  beginning 
column  number  and  then  an  ending  column  number.  The  remainder  of  the 
record  would  contain  the  grid  matrix  numbers  corresponding  to  elements 
defined  from  (row,  beginning  column)  to  (row,  ending  column) .  An 
example  of  this  file  structure  is  provided  on  the  following  page. 

Before  beginning  the  grid-data  input,  the  file  must  contain  two 
header  records.  These  records  can  be  blank  but  typically  contain 


GWPATH  Version  3.0,  1987 
OUTPUT  CONTROL  INFORMATION 


Printable  File  of  Results  [EXAMPLE.DAT  ] 

Press  ENTER  If  Option  Not  Desired 

Plot  Results  With  HP7475A/IBM7372  Color  Plotter,  (Y)es  or  (N)o? 


Figure  7d.  Screen  4:  Output  control  menu 
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Illinois  State  Water  Survey 


2204  Griffith  Drive 
Champaign,  Illinois  61820-7495 


ORDER  FORM 


To  receive  a  S?",  1.2  Mb  diskette  containing  GWPATH3.0,  IBM 
device  drivers,  and  example  data  sets,  please  send  check  (or 
money  order)  for  $125.00  payable  to  University  of  Illinois,  c/o 
Illinois  State  Water  Survey  to  John  M.  Shafer  at  the  above 
address. 


Name  _ 

Organization  _ 

Address  _ 

City  _  State _ Zip 

Date  Phone 


(Please  allow  2  to  3  weeks  for  delivery) 


information  regarding  the  flow  domain.  The  following  FORTRAN  source  code 
is  representative  of  how  the  files  are  read  by  GWPATH. 

OPEN(IUNIT,FILE=' filename. ext' ) 

READ (lUNIT, 100, END=300)  NC ,  NR,  NULL,  XMIN,  XMAX,  YMIN,  YMAX, 

ZMIN,  ZMAX 

100  FORMAT(3I5,6F10.3) 

READ (lUNIT, 100, END=300)  NULL,  NULL,  NULL,  ZMIN,  ZMAX 
200  READ (lUNIT, 100, END=300)  J,  II,  12,  (Z(J,I),  I  =  11,12) 

GO  TO  200 
300  CLOSE  (lUNIT) 


An 

example  ASCII 

data  file 

with  the 

proper  structure  for 

GWPATH 

input 

;  is 

shown 

using  the 

hydraulic 

conductivity  data 

listed  below, 

5  1 

5  1 

0 

0.000 

8000.000 

0.000 

8000.000 

10.000 

20.000 

0 

0 

0 

10.000 

20.000 

1 

1 

6 

11.258 

1 1 . 686 

1  1  .887 

12.025 

12.221 

12.513 

1 

7 

1  2 

1 2 . 792 

13.010 

13.168 

13.260 

13.280 

13.227 

1 

1  3 

1  8 

13.104 

12.923 

12.700 

12.450 

12.184 

11.913 

1 

1  9 

24 

11.645 

11.386 

11.142 

10.917 

10.713 

10.534 

1 

25 

30 

10.379 

10.251 

10.149 

10.074 

10.024 

10.000 

1 

31 

36 

10.000 

10.024 

10.070 

10.136 

10.221 

10.323 

1 

37 

42 

10.440 

10.569 

10.708 

10.853 

11.000 

11.148 

1 

43 

48 

11.291 

11.425 

11.548 

11.653 

11.737 

1 1 . 796 

1 

49 

5  1 

11  .826 

11.842 

11.858 

2 

1 

6 

11.148 

11.580 

11.954 

12.176 

12.405 

1 2 .666 

2 

7 

1  2 

12.951 

13.189 

13.370 

13.494 

13.553 

13.544 

2 

1  3 

1  8 

13.507 

13.340 

13.087 

1 2 . 868 

12.628 

12.376 

2 

1  9 

24 

12.117 

1 1 .866 

11.626 

11.403 

11.200 

11.019 

2 

25 

30 

10.861 

10.727 

10.619 

10.535 

10.477 

10.442 

2 

3  1 

36 

10.430 

10.441 

10.472 

10.522 

10.591 

10.674 

2 

37 

42 

10.771 

10.880 

1 0 .997 

11.119 

11.244 

11.368 

2 

43 

48 

11.488 

11.601 

11.702 

1 1 . 788 

11.856 

11.902 

2 

49 

5  1 

11.925 

11.940 

11.954 

5  1 

1 

6 

16.254 

16.224 

16.200 

16.183 

16.170 

16.141 

5  1 

7 

1  2 

16.105 

16.081 

16.079 

16.102 

16.139 

16.180 

5  1 

1  3 

1  8 

16.224 

16.270 

16.318 

16.367 

16.416 

16.465 

5  1 

1  9 

24 

16.513 

16.561 

1  6 .608 

16.654 

1 6 .698 

16.741 

5  1 

25 

30 

1 6 . 783 

16.824 

1 6 .863 

16.901 

16.938 

16.973 

5  1 

3  1 

36 

17.007 

1 7.038 

1 7. 068 

17.096 

17.120 

17.142 

5  1 

37 

42 

17.159 

17.172 

17.179 

17.180 

17.172 

17.155 

5  1 

43 

48 

17.128 

17.087 

1 7.033 

1 6 . 963 

16.879 

16.781 

5  1 

49 

5  1 

1 6 .669 

16.543 

16.404 

EOF 


19 


GWPATH  also  allows  an  optional  ASCII  data  file  containing  the 
ground-water  flow  domain  source/sink  location  information.  This  file 
must  also  be  prepared  for  list- directed  input.  The  first  record  of  the 
source/sink  file  must  contain  the  number  of  sources  followed  by  the 
number  of  sinks.  The  real  x-y  coordinate  locations  of  all  sources  are 
listed  first  followed  by  the  real  x-y  coordinate  locations  of  all  known 
sinks.  The  following  example  shows  how  the  ASCII  source/sink  data  file 
must  be  organized. 

Record  No . 

1 
2 
3 


NSOUR+1 

NSOUR-l-2 

NSOUR-l-3 


NSOUR+NSINK-t-1 

EOF 


Output  Specifications 
The  user  is  also 

the  results  of  the  pathline  analysis.  If  hardcopy  output  is  not 
required  for  any  particular  analysis  session,  pressing  <Enter>  before 
any  other  response  will  eliminate  hardcopy  output.  If  a  filename  is 
provided,  GWPATH  writes  the  results  (i.e.,  the  pathline(s)  trajectory 
and  travel  time)  to  this  file  which  can  be  printed  after  the  program 
ends . 


Value 

NSOUR  NSINK 

SOURCE (l,x)  SOURCE (l,y) 
SOURCE ( 2, x)  SOURCE ( 2, y) 


SOURCE (NSOUR, x)  SOURCE (NSOUR, y) 
SINK(l,x)  SINK(l,y) 

SINK(2,x)  SINK(2,y) 


SINK(NSINK,x)  SINK(NSINK,y) 


prompted  to  name  an  ASCII  disk  file  for  storing 


The  final  GWPATH  control  information  input  requirement  is  a  "Yes" 
or  "No"  response  to  whether  or  not  a  hardcopy  plot  (i.e.,  HP-7475A  or 
IBM-7372  six-pen  plotter  output)  is  desired.  If  the  response  to  this 
question  is  "Yes",  the  plot  will  mirror  the  screen  graphics  until 
program  execution  is  halted. 


The  following  pen  stall  numbers  correspond  to  the  indicated 
components  of  a  hardcopy  plot. 


Stall  No. 

1 

2 

3 

4 

5 

6 


Component  of  Plot 

Flow  domain  border, 
legend  box 
Text,  marker  symbol 
outline 

Flow  domain  grid 
Sources ,  forward 
tracked  paths 
Sinks,  reverse 
tracked  paths 
Capture  zone 


Recommended  Pen 

0.7mm  Black 

0 . 3mm  Black 

0 . 3mm  Yellow 
0 . 3mm  Blue 

0 . 3mm  Red 

0 . 7mm  Blue 
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GWPATH  Execution  Control 


GWPATH  execution  control  enables  the  user  to  interactively  specify 
several  options  for  running  the  program,  especially  the  manner  in  which 
pathline  starting  locations  are  determined.  At  the  beginning  of  the 
pathline  calculation  procedure,  a  graphical  representation  of  the  flow 
domain  is  displayed  on  the  screen  followed  by  a  command  line  displayed 
at  the  bottom.  Figure  8  shows  the  display  screen  arrangement  of  the 
flow  domain,  legend,  time  -  stepping  control  information,  and  the  command 
line.  The  user  must  respond  via  keyboard  input  to  the  various  options 
available  for  program  execution. 

The  first  option  to  choose  is  whether  pathline  analysis  or  capture 
zone  analysis  is  to  be  performed.  Pathline  analysis  allows  for  forward 
or  reverse  tracking,  multiple  starting  locations,  and  display  of 
pathlines.  Capture  zone  analysis  automatically  invokes  reverse  tracking 
and  suspends  display  of  individual  pathlines.  However,  capture  zone 
analysis  enables  the  feature  that  connects  the  end  points  of  all 
pathlines  and  will  cross-hatch  the  interior  of  the  resulting  polygon. 

If  pathline  analysis  is  the  option  selected,  the  user  is  asked  if 
forward  pathline  tracking  or  reverse  pathline  tracking  is  to  be 
performed.  The  next  three  questions  relate  to  the  time  specifications 
for  the  analysis.  The  first  question  asks  for  the  maximum  allowable 
travel  time,  in  days.  After  the  maximum  travel  time  is  entered,  the 
user  is  asked  for  the  minimum  number  of  moves  per  grid  cell  and  the 
absolute  maximum  time  step,  in  days.  The  previous  discussion  on  time¬ 
stepping  explains  how  these  parameters  affect  the  calculation  of 
pathlines . 

Once  the  time  -  stepping  specifications  have  been  entered,  the  user 
is  asked  to  choose  the  mode  of  pathline  starting  location  data  entry. 

The  five  different  ways  to  enter  pathline  starting  locations  are: 

1.  Manual  cursor  entry  (pathline  analysis  only) 

2.  Manual  keyboard  entry  (pathline  analysis  only) 

3.  Circular  starting  configuration  from  cursor  entry 

4.  Circular  starting  configuration  from  keyboard  entry 

5.  ASCII  file  of  starting  locations  (pathline  analysis  only) 

If  cursor  entry  is  selected,  a  cross  hair  appears  on  the  screen  in  the 
middle  of  the  flow  domain.  The  arrow  keys  are  used  to  move  the  cursor 
to  the  desired  pathline  starting  location.  The  <Insert>  key  toggles 

between  large  cursor  increments  and  small  cursor  increments .  At  any 
time  the  real  coordinate  location  can  be  determined  by  pressing  <Enter>. 
To  select  any  particular  location  to  start  a  pathline  press  <S>. 

Multiple  pathlines  can  be  calculated  by  repeating  this  process  a  maximum 
of  299  times.  To  end  selecting  pathline  starting  locations  and  begin 
their  calculation  the  user  must  press  <Esc>. 

To  select  precise  starting  locations,  the  user  should  use  the 
keyboard  input  option.  With  this  option  the  user  types  the  exact 
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Figure  8 .  Example  display  screen  for  GWPATH 


(integer)  x-y  coordinates  of  pathline  starting  locations.  In  either 
case,  once  all  starting  locations  have  been  entered  the  user  is  asked  if 
the  locations  should  be  saved  in  an  ASCII  file.  If  the  answer  is  "Yes" 
the  user  must  enter  a  filename  for  the  locations. 

Instead  of  entering  each  pathline  starting  location  individually, 
pathline  starting  locations  can  be  automatically  calculated  around  a 
circle  defined  by  the  user.  This  option  requires  the  user  to  enter  the 
x-y  coordinates  of  the  center  of  a  circle  and  then  provide  a  radius  from 
the  center.  GWPATH  will  then  automatically  calculate  pathline  starting 
locations  at  equal  arc  lengths  around  the  circle.  As  is  the  case  for 
manual  starting  location  entry,  the  center  of  the  circle  can  be  input  by 
using  the  cursor  or  by  direct  keyboard  entry. 
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In  addition  to  entering  pathline  starting  locations  either 
manually  or  with  a  circular  configuration,  GWPATH  allows  for  starting 
location  input  from  an  ASCII  file.  If  the  ASCII  disk  file  option  is 
selected  for  input  of  pathline  starting  locations,  the  file  must  be 
formatted  for  list-directed  input  according  to  the  following  example. 


Record  No. 

1 

2 


Value 

BECIN(l,x)  BECIN(l,y) 
BECIN(2,x)  BECIN(2,y) 


NPATH  BEGIN (NPATH,x)  BEGIN (NPATH , y) 

EOF 


GWPATH  Dimensionality 


GWPATH  dimension  limitations  are  as  follows: 


•  Maximum  number 

•  Maximum  number 

•  Maximum  number 

•  Maximum  number 

•  Maximum  number 

•  Maximum  number 


of  columns  =  99 
of  rows  =  99 
of  sources  =  20 
of  sinks  =  20 

of  pathlines  per  execution  =  300 
of  discrete  points  per  pathline  =  1000 
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OUTPUT 


By  entering  a  filename  at  the  "Printable  File  of  Results  [  ]" 
prompt,  an  ASCII  file  of  the  pathline  trajectories  and  travel  times  will 
be  created.  The  file  will  automatically  be  saved  in  the  current 
directory  when  execution  of  GWPATH  is  terminated.  The  file  can  then  be 
printed  in  the  same  manner  as  any  ASCII  file.  If  <Enter>  is  pressed  at 
the  prompt,  the  results  of  the  current  session  will  not  be  saved. 


following 

Point 

listing  is  an  example  of  GWPATH  file  output 

DEMONSTRATION  OF  OUTPUT  TO  FILE 

Pathline  No.  =  1 

Travel  Time  (days)  =  500.00 

No.  of  Points  =  51 

X  (  feet  )  Y  (  feet  )  Time  (davs) 

1 

1943.00 

1814.00 

0.00 

2 

1946.24 

1823.08 

10.00 

3 

1949.45 

1832.66 

20.00 

4 

1952.53 

1842.73 

30.00 

5 

1955.35 

1853.23 

40.00 

11 

1961.26 

1928.26 

100.00 

12 

1960.82 

1943.37 

110.00 

13 

1960.66 

1959.35 

120.00 

14 

1960.77 

1976.26 

130.00 

15 

1961.17 

1994.15 

140.00 

25 

1979.11 

2161.50 

240.00 

26 

1981.47 

2176.13 

250.00 

27 

1983.87 

2190.54 

260.00 

28 

1986.29 

2204.74 

270.00 

29 

1988.75 

2218.73 

280.00 

30 

1991.23 

2232.52 

290.00 

45 

2030.26 

2419.51 

440.00 

46 

2032.92 

2431.10 

450.00 

47 

2035.60 

2442.60 

460.00 

48 

2038.30 

2454.01 

470.00 

49 

2041.01 

2465.33 

480.00 

50 

2043.74 

2476.57 

490.00 

51 

2046.47 

2487.74 

500.00 

Ik 


GWPATH  VALIDATION 


GWPATH  was  validated  by  forward  tracking  a  pathline  and  comparing 
the  estimated  trajectory  with  a  known  calculated  trajectory.  The 
reverse  tracking  option  was  subsequently  validated  by  reverse  tracking 
the  known  pathline.  Finally,  the  effective  porosity  was  divided  by  2.0, 
doubling  the  velocity.  Therefore,  it  should  take  one -half  the  original 
travel  time  for  the  same  trajectory,  all  other  factors  remaining  the 
same . 


To  accomplish  the  validation  a  regular  8000  ft  by  8000  ft 
homogeneous,  isotropic  flow  domain  was  constructed  with  a  hydraulic  head 
distribution  as  shown  in  figure  9.  The  validation  flow  domain  is 
symmetric  about  the  north-south  centerline;  i.e.,  x  =  4000  ft. 

Consequently,  pathlines  along  the  centerline  are  linear.  This 
characteristic  is  used  to  validate  GWPATH.  The  following  parameters 
describe  the  validation  flow  domain  and  the  imposed  aquifer  properties. 

•No.  of  column  nodes  =  51 
•No.  of  row  nodes  =  51 

•  Delta-X  =  160  ft 

•  Delta-Y  =  160  ft 

•  Kjj  =  50  ft/day 

•  Ky  =  50  ft/day 

•  ne  =  0.5 

The  following  time  -  stepping  constraints  were  used  for  the  code 
validation. 

•  Maximum  allowable  travel  time  =  10,000  days 

•  Minimum  number  of  moves  per  grid  cell  =  1 

•  Maximum  allowable  time  step  =  25  days 
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Figure  9.  Hydraulic  head  (feet  above  MSL)  for  GWPATH  validation 
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The  forward  tracking  validation  pathline  computed  by  GWPATH  begins 
at  location  (4000  ft,  1120  ft).  The  terminating  coordinate  (i.e., 
s[10,000  days])  equals  (4000  ft,  7352.87  ft).  The  pathline  trajectory 
is  presented  in  the  following  listing  and  is  shown  in  figure  10. 

FORWARD  TRACKING  VALIDATION 


Pathline  No.  =  1 
Travel  Time  (days)  =  10000.00 
No.  of  Points  =  401 


lint 

X  (  feet  ) 

Y  (  feet  ) 

Time  (dav: 

1 

4000.00 

1120.00 

0.00 

2 

4000.00 

1136.09 

25.00 

3 

4000.00 

1152.17 

50.00 

4 

4000.00 

1168.25 

75.00 

5 

4000.00 

1184.33 

100.00 

6 

4000.00 

1200.41 

125.00 

7 

4000.00 

1216.49 

150.00 

8 

4000.00 

1232.56 

175.00 

9 

4000.00 

1248.64 

200.00 

10 

4000.00 

1264.71 

225.00 

11 

4000.00 

1280.78 

250.00 

12 

4000.00 

1296.85 

275.00 

13 

4000.00 

1312.92 

300.00 

14 

4000.00 

1328.98 

325.00 

15 

4000.00 

1345.05 

350.00 

16 

4000.00 

1361.11 

375.00 

17 

4000.00 

1377.17 

400.00 

18 

4000.00 

1393.23 

425.00 

19 

4000.00 

1409.29 

450.00 

385 

4000.00 

7.110.79 

9600.00 

386 

4000.00 

7125.93 

9625.00 

387 

4000.00 

7141.07 

9650.00 

388 

4000.00 

7156.21 

9675.00 

389 

4000.00 

7171.35 

9700.00 

390 

4000.00 

7186.48 

9725.00 

391 

4000.00 

7201.61 

9750.00 

392 

4000.00 

7216.75 

9775.00 

393 

4000.00 

7231.88 

9800.00 

394 

4000.00 

7247.01 

9825.00 

395 

4000.00 

7262.13 

9850.00 

396 

4000.00 

7277.26 

9875.00 

397 

4000.00 

7292.38 

9900.00 

398 

4000.00 

7307.51 

9925.00 

399 

4000.00 

7322.63 

9950.00 

400 

4000.00 

7337.75 

9975.00 

401 

4000.00 

7352.87 

10000.00 
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Figure  10.  Forward  tracked  GWPATH  validation  pathline 


To  validate  the  code,  this  numerically  determined  pathline  was 
analytically  derived.  The  average  hydraulic  gradient  over  the  length  of 
the  pathline  was  calculated  by  taking  the  mean  of  the  hydraulic  gradient 
in  the  north-south  direction  between  each  pair  of  nodes  in  the  flow 
domain  grid  between  the  starting  point  and  the  ending  point  of  the 
pathline.  The  resulting  average  gradient  is 

^  =  0.0062339 


With  equation  3,  the  average  pore  ground-water  velocity 
trajectory  of  the  pathline  can  be  calculated  as 


X  0.0062339  =  0.62339  ft/day 


along 


the 


Multiplying  v^  by  10,000  days  results  in  a  linear  travel  distance  of 
6,233.9  ft.  Adding  the  travel  distance  to  the  starting  location  gives 
the  pathline  ending  location  (4000  ft,  7353.9  ft).  Comparison  of  this 
value  with  the  ending  location  of  the  validation  pathline  calculated  by 
GWPATH  shows  that  within  roundoff  error  the  code  correctly  estimates 
ground-water  pathlines  in  the  forward  tracking  mode.  The  error  between 
the  GWPATH  numerically  calculated  pathline  and  the  same  pathline 
analytically  determined  is  approximately  0.01%. 


After  the  forward  pathline  tracking  methodology  implemented  in 
GWPATH  had  been  validated,  the  reverse  tracking  was  easily  validated  by 
simply  reverse- tracking  the  forward  tracked  pathline.  The  same  time 
stepping  constraints  used  for  the  forward  pathline  validation  were  used 
for  the  reverse  pathline  validation.  The  pathline  starting  location  for 
reverse  pathline  validation  was  the  ending  location  of  the  forward 
tracked  pathline;  i.e.,  (4000  ft,  7352.87  ft).  The  total  allowable 
travel  time  was  10,000  days.  The  following  pathline  trajectory  is  the 
result  of  the  reverse  tracking  validation.  The  pathline  for  reverse 
calculation  validation  is  presented  in  figure  11. 
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REVERSE  TRACKING  VALIDATION 


Pathline  No.  =  1 
Travel  Time  (days)  =  10000.00 
No.  of  Points  =  401 


X  (  feet  ) 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 


Y  (  feet  ) 

7352.87 

7337.75 
7322.63 

7307.51 
7292.39 
7277.26 
7262.14 
7247.01 

7231.88 

7216.75 
7201.62 
7186.48 
7171.35 
7156.21 
7141.08 
7125.94 
7110.80 
7095.65 

7080.51 


Time  (days) 

0.00 

25.00 

50.00 

75.00 

100.00 

125.00 

150.00 

175.00 

200.00 

225.00 

250.00 

275.00 

300.00 

325.00 

350.00 

375.00 

400.00 

425.00 

450.00 


4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 

4000.00 


1377.18 

1361.11 

1345.05 

1328.99 

1312.92 

1296.85 

1280.78 

1264.71 

1248.64 

1232.57 

1216.49 

1200.41 

1184.33 

1168.25 

1152.17 

1136.09 

1120.00 


9600.00 

9625.00 

9650.00 

9675.00 

9700.00 

9725.00 

9750.00 

9775.00 

9800.00 

9825.00 

9850.00 

9875.00 

9900.00 

9925.00 

9950.00 

9975.00 

10000.00 
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Figure  11.  Reverse  tracked  GWPATH  validation  pathline 


The  results  of  the  reverse  tracking  exercise  clearly  validate  the 
reverse  pathline  analysis  capability  of  GWPATH.  The  pathline  ending 
location  ( in  a  backward  sense)  after  10,000  days  is  exactly  the  starting 
location  for  the  same  pathline  forward  tracked;  i.e.,  (4000  ft,  1120 
ft). 

The  final  code  validation  exercise  was  to  decrease  the  effective 
porosity  by  a  factor  of  2,  thereby  doubling  the  average  pore  ground- 
water  velocity.  Instead  of  0.50,  ng  equalled  0.25  for  this  case.  The 
time  stepping  constraints  remained  the  same  as  for  the  previous 
validation  exercises.  However,  the  total  allowable  travel  time  was  also 
reduced  by  a  factor  of  2  to  reconstruct  the  original  pathline.  The  same 
forward  tracking  pathline  starting  location  was  used;  i.e.,  (4000  ft, 
1120  ft) .  The  anticipated  results  were  that  the  length  of  the  pathline 
calculated  by  GWPATH  should  be  identical  to  the  previous  pathlines 
because  the  mean  velocity  was  doubled  while  the  travel  time  was  halved. 
The  resulting  pathline  trajectory  is  as  follows. 
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REDUCED  EFFECTIVE  POROSITY 


Pathline  No.  =  1 
Travel  Time  (days)  =  5000.00 
No.  of  Points  =  201 


oint 

X  (  feet  ) 

Y  (  feet  ) 

Time  (davs 

1 

4000.00 

1120.00 

0.00 

2 

4000.00 

1152.17 

25.00 

3 

4000.00 

1184.33 

50.00 

4 

4000.00 

1216.49 

75.00 

5 

4000.00 

1248.64 

100.00 

6 

4000.00 

1280.78 

125.00 

7 

4000.00 

1312.92 

150.00 

8 

4000.00 

1345.05 

175.00 

9 

4000.00 

1377.17 

200.00 

10 

4000.00 

1409.29 

225.00 

191 

4000.00 

7050.21 

4750.00 

192 

4000.00 

7080.50 

4775.00 

193 

4000.00 

7110.79 

4800.00 

194 

4000.00 

7141.07 

4825.00 

195 

4000.00 

7171.35 

4850.00 

196 

4000.00 

7201.61 

4875.00 

197 

4000.00 

7231.88 

4900.00 

198 

4000.00 

7262.13 

4925.00 

199 

4000.00 

7292.38 

4950.00 

200 

4000.00 

7322.63 

4975.00 

201 

4000.00 

7352.87 

5000.00 

As  expected,  the  pathline  is  identical  to  the  previous  pathline,  only 
the  travel  time  is  one -half  the  original  10,000  days. 


30 


GWPATH  DEMONSTRATION 


Several  example  applications  of  GWPATH  are  presented  in  order  to 
demonstrate  its  use.  For  demonstration  purposes,  a  complex  ground- 
water  flow  domain  was  developed  with  the  following  specifications. 


•No.  of  column  nodes  =  51 
•No.  of  row  nodes  =  51 

•  Delta-X  =  160  ft 

•  Delta-Y  =  160  ft 

•  Inhomogeneous,  anisotropic  aquifer  conditions 

•  Constant  head  boundaries 

•  Steady  state  hydraulic  head 

•  Five  hydraulic  sources;  five  hydraulic  sinks 


Figure  12  is  a  contour  plot  of  the  steady  state  hydraulic  head 
distribution  of  the  demonstration  ground-water  flow  domain.  Figure  13 
shows  the  spatial  variability  of  x-direction  hydraulic  conductivity 


Figure  12.  Hydraulic  head  (feet  above  MSL)  for  GWPATH  flow  path 
analysis  demonstration 
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which  ranges  from  10  ft/day  to  20  ft/day.  In  similar  fashion,  figure  14 
shows  the  y-direction  spatially  varying  hydraulic  conductivity  which 
ranges  from  10  ft/day  to  200  ft/day.  Finally,  the  effective  porosity  is 
assumed  to  range  from  a  low  of  0.15  to  a  high  of  0.25  over  the  flow 
domain  (figure  15). 


Figure  13.  X-direction  hydraulic  conductivity  (ft/day)  for  GWPATH 
flow  path  analysis  demonstration 
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Figure  14.  Y-direction  hydraulic  conductivity  (ft/day)  for  GWPATH 
flow  path  analysis  demonstration 
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Figure  15.  Effective  porosity  (dimensionless)  for  GWPATH  flow  path 
analysis  demonstration 
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The  five  hydraulic  sources  and  their  steady  state  source  strengths 


X- Coordinate 

1760  ft 
3680  ft 
1120  ft 
1440  ft 
2880  ft 


Y-Coordinate 

1120  ft 
1760  ft 
2240  ft 
3200  ft 
3200  ft 


Injection  Rate 

200  gpm 
400  gpm 
600  gpm 
500  gpm 
500  gpm 


The  five  hydraulic  sinks  and  their  steady  state  source  strengths  are: 


X-Coordinate  Y-Coordinate  Pumping  Rate 


3040 

ft 

4640  ft 

400 

gpm 

4320 

ft 

4320  ft 

400 

gpm 

5440 

ft 

3680  ft 

400 

gpm 

4640 

ft 

6560  ft 

500 

gpm 

5760 

ft 

5600  ft 

600 

gpm 

Figure  16  shows  the  orientation  of  sources  and  sinks  throughout  the  flow 
domain . 
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Figure  16.  Source/sink  locations  for  GWPATH  flow  path  analysis 
demonstration 
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Example  No.  1:  Forward  Tracking 

This  example  demonstrates  the  forward  pathline  tracking  feature  of 
GWPATH.  Both  the  manual  selection  and  circular  configuration  options 
for  pathline  starting  locations  were  used  to  perform  the  pathline 
analysis  shown  in  figure  17.  The  pathlines  beginning  in  the  upper 
righthand  (i.e.,  northeast)  boundary  of  the  flow  domain  were 
individually  selected  using  the  cursor  input  option.  The  starting 
locations  of  the  forward  tracked  pathlines  emanating  from  the  source 
located  at  (2880  ft,  3200  ft)  were  calculated  using  the  circular 
starting  location  configuration  option.  Equal  arc  distances  around  an 
imaginary  circle  centered  at  (2880  ft,  3200  ft)  with  a  100  ft  radius 
were  used  to  automatically  compute  20  pathline  starting  locations.  The 
following  time  -  stepping  constraints  were  used  for  this  example. 

•  Maximum  travel  time  =  5000  days 

•  Minimum  moves  per  cell  =  2 

•  Maximum  allowable  time  step  =  25  days 

Comparison  of  the  pathlines  calculated  by  GWPATH  (figure  17)  with  the 
hydraulic  head  distribution  (figure  12)  confirms  the  sensitivity  of  the 
methodology  to  complex  flow  domains. 


0  0000 
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Figure  17.  GWPATH  forward  pathline  tracking  demonstration 
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Example  No.  2:  Reverse  Tracking 

To  demonstrate  the  reverse  pathline  tracking  option  of  GWPATH ,  a 
circular  pathline  starting  configuration  around  a  hydraulic  sink  was 
employed.  In  a  backward  sense,  however,  the  input  pathline  starting 
locations  are,  in  reality,  the  end  points  of  the  pathlines  in  the 
direction  of  flow.  GWPATH  simply  tracks  these  pathlines  up-gradient 
(i.e.,  backward)  as  opposed  to  forward,  or  down- gradient ,  through  the 
flow  domain.  Gonsequently ,  instead  of  being  convergent  toward  a 
hydraulic  sink,  the  reverse  pathline  tracking  flow  paths  are  divergent 
away  from  a  sink. 

The  same  flow  domain  used  for  Example  No.  1  (figures  12  -  16)  was 
used  for  the  reverse  tracking  example  application.  Twenty  equally 
spaced  pathline  starting  locations  were  automatically  calculated  by 
GWPATH  along  an  imaginary  circle  of  160  ft  radius  around  the  sink 
located  at  (4640  ft,  6560  ft).  The  following  time  -  stepping  constraints 
were  imposed  on  the  analysis. 

•  Maximum  travel  time  =  10,000  days 

•  Minimum  moves  per  cell  =  2 

•  Maximum  allowable  time  step  =  25  days 

The  results  of  the  reverse  tracking  demonstration  are  presented  in 
figure  18.  The  ability  of  GWPATH  to  adequately  evaluate  complex  flow 
domains  is  well  represented  in  this  figure.  The  pathlines  calculated  by 
GWPATH  and  displayed  in  figure  18  are  not  necessarily  obvious  upon 
visual  inspection  of  the  hydraulic  head  contour  map  presented  in  figure 
12.  However,  the  reverse  pathline  analysis  clearly  shows  the  hydraulic 
relationship  of  the  sink  at  (4640  ft,  6560  ft)  to  the  hydraulic  sources 
in  the  flow  domain  and  to  the  other  sinks.  Further,  given  that  by 
definition  no  flow  can  be  orthogonal  to  a  pathline,  the  contribution 
from  various  sources  of  flow  can  be  estimated  for  any  given  sink  from 
the  distribution  of  reverse  pathlines  emanating  from  the  sink.  If,  at  a 
hydraulic  sink,  reverse  tracked  pathlines  are  started  along  even  arc 
distances  around  the  sink  (as  was  the  case  for  this  example)  ,  there  is 
equal  flow  contribution  in  each  flow  tube  represented  by  any  adjacent 
pair  of  pathlines. 

Example  No.  3:  Capture  Zone  Analysis 

In  addition  to  forward  and  reverse  pathline  tracking,  an  extension 
of  reverse  tracking  is  included  in  GWPATH  which  allows  the  calculation 
of  time-related  capture  zones  in  ground-water  flow  domains.  For  a 
complete  discussion  of  reverse  pathline  calculation  of  time-related 
capture  zones  using  GWPATH,  refer  to  Shafer  (1987). 

A  t-year  capture  zone  in  a  three-dimensional  flow  regime  is  the 
bounded  volume  described  by  the  set  of  all  flow  paths  that  intercept  a 
hydraulic  sink  in  t  years.  For  the  case  of  two-dimensional  horizontal 
flow  the  vertical  flow  component  is  ignored  and  the  capture  zone  becomes 
the  planar  area  described  by  the  set  of  all  horizontal  pathlines  with  t- 
year  travel  times  to  the  sink. 
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Figure  18.  GWPATH  reverse  pathline  tracking  demonstration 


GWPATH  estimates  the  time -related  capture  zone  by  reverse  tracking 
a  large  number  of  pathlines  (typically  between  100  and  300)  all  with  the 
same  travel  time.  GWPATH  connects  the  resulting  end  points  to  form  a 
polygon  around  the  sink  which  represents  the  t-year  capture  zone. 
GWPATH  then  fills  the  polygon  with  a  hatch  pattern  to  display  the 
enclosed  area. 


The  flow  domain  represented  by  the  steady  state  hydraulic  head 
distribution  shown  in  figure  19  was  used  by  GWPATH  to  calculate  a  time- 
related  capture  zone  for  demonstration.  The  flow  domain  was  assumed 
to  be  inhomogeneous  and  isotropic  with  block  hydraulic  conductivities 
shown  in  figure  20.  Effective  porosity  was  assumed  to  be  a  constant 
0.25.  Four  steady  strength  hydraulic  sinks  were  located  in  the  flow 
domain  at  the  following  coordinate  locations. 
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The  10-year  (i.e.,  3652  day)  capture  zone  was  calculated  for  the 
hydraulic  sink  located  at  (3840  ft,  4640  ft).  The  following  time 
stepping  constraints  were  used  for  the  analysis. 
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Figure  19.  Hydraulic  head  (feet  above  MSL)  for  GWPATH  capture  zone 
analysis  demonstration 


•  Maximum  travel  time  =  3652  days 

•  Minimum  moves  per  cell  =  2 

•  Maximum  allowable  time  step  =  20  days 

For  this  demonstration,  300  reverse  tracked  pathlines  were  used  to 
delimit  the  capture  zone.  Figure  21  shows  the  resulting  10-year  capture 
zone  calculated  by  GWPATH  around  the  sink  at  (3840  ft,  4640  ft). 
Although  not  performed  for  this  example,  once  the  capture  zone  has  been 
calculated,  the  flow  path  analysis  mode  can  be  re-entered  and  a  limited 
number  of  reverse  pathlines  (e.g.,  20)  overplotted  on  the  capture  zone 

for  verification  purposes  or  to  estimate  flow  contributions  from 
specific  points  in  the  flow  domain. 
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Figure  20.  Block  hydraulic  conductivities  (ft/day)  for  GWPATH 
capture  zone  analysis  demonstration 
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Figure  21.  10-year  capture  zone  for  hydraulic  sink  located 
at  (3840  ft,  4640  ft)  for  GWPATH  demonstration 
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APPENDIX 

Software  Installation  Instructions 
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Installing  GWPATH 


1.  Establish  directory  on  hard  disk  for  GWPATH. 

mkdir  c:\gwpath 

2.  Change  directory. 

chdir  c:\gwpath 

3.  Place  diskette  in  drive  A. 

4.  Copy  diskette  to  hard  disk. 

copy  a:\*.'^  /v 


Running  GWPATH 


To  run  GWPATH  type  GWPATH 

NOTE:  All  data  files  must  reside  in  the  current  directory  from  which  the 
call  to  run  GWPATH  is  made. 
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